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Abstract 

An asymptotic theory is developed to generate equations that model the global behaviour of electromagnetic waves 
in periodic photonic structures when the wavelength is not necessarily long relative to the periodic cell dimensions; 
potentially highly-oscillatory short-scale detail is encapsulated through integrated quantities. Our approach is 
based on the method of high-frequency homogenisation (HFH), introduced in the context of scalar out-of-plane 
elastic waves in pO) , here extended to the full Maxwell system for electromagnetic waves in three dimensions. 
In doing so, a vector treatment of Maxwell’s equations is necessary, of which scalar transverse electric (TE) and 
transverse magnetic (TM) polarised cases are a strict subset. The resulting procedure yields effective dynamic 
continuum equations, valid even at high frequencies, that govern a scalar function providing long-scale modulation 
of short-scale Bloch eigensolutions inside the structure. The form of the effective equation changes in the case of 
non-trivially repeated eigenvalues, leading to degeneracies that are linked to the appearance of Dirac-like points, a 
case that we explore using the asymptotic method. We examine the low-frequency long-wave limit in the case of 
dielectric structures, and find the expected result of classical homogenisation theory is captured, in which the usual 
governing equation is recovered, but now with the permittivity replaced by an effective tensor. 

The theory we develop is then applied to two topical examples, the first being the case of aligned dielectric cylin¬ 
ders, which has great importance in the modelling of photonic crystal fibres. Advantage can be taken of the 
axial invariance to reduce the three-dimensional example such that it permits manageable computations on a two- 
dimensional domain. Results of the asymptotic theory are verified against these numerical simulations by compar¬ 
ing photonic band diagrams and evanescent decay rates for localised modes. We then consider the propagation of 
waves in a structured metafilm, here chosen to be a planar array of dielectric spheres. These waves can be highly 
localised in the plane, decaying exponentially in directions normal to the array, and at certain frequencies strongly 
directional dynamic anisotropy is observed. Computationally this is a challenging three-dimensional calculation, 
which we do here, and then demonstrate that the asymptotic theory captures the effect, giving highly accurate 
qualitative and quantitative comparisons as well as providing interpretation for the underlying change from elliptic 
to hyperbolic behaviour. 


1 Introduction 

The practical implementation of photonic crystal fibres (PCFs) in communications, lasers and other areas of en¬ 
gineering has resulted in much interest being ascribed to the dispersive properties of periodic photonic struc- 
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tures 134 in particular those exhibiting photonic band gaps in one or more directions. Further to this, the rich 
band structures possessed by photonic crystals (PCs) give rise to a host of novel physical phenomena including, but 
not limited to, dynamic anisotropy fTS) , ultrarefraction Q and light confinement |j^ that can also impact upon the 
design of PCFs | [53] . From a numerical point of view, the implementation of Maxwell’s equations in such struc¬ 
tures can pose considerable difficulty, particularly when the wavelength is of a similar scale to that of the repeating 
microcell, of which there may be many hundreds or more. In such cases, multiple scattering and interference play 
a central role. Typical numerical schemes such as finite element methods | |62| are highly effective, but the compu¬ 
tational cost in tackling the coupled vector system increases rapidly with the number of nodes, which is intimately 
connected with the dimensionality, scale and operating wavelength of the system. Optimisation and design of a PC 
or PCF may require multiple simulations across changes in geometry, material parameters and microstructure, so 
it is attractive to explore complementary effective medium, or homogenisation, schemes 00 that attempt to treat 
the medium in some averaged sense and use these instead of, or alongside, direct numerical simulation. 

Asymptotic homogenisation dates back to the French, Italian and Russian schools, and was developed for studying 
Partial Differential Equations (PDEs) with rapidly oscillating coefficients. The approach of introducing indepen¬ 
dent fast and slow variables applies naturally to structured media, leading to ‘effective medium’ models on which 
there is extensive literature ||^|9 The effective properties attained using such methodologies are valid 

within the low-frequency (quasi-static) regime, which is the relevant domain for the rapidly developing theory of 
sub-wavelength metamaterials, in which a typical application involves tuning these effective properties to desired 
values as prescribed by transformation optics (see, for example | |^[^ ). Unfortunately these homogenisation the¬ 
ories often perform poorly near band edges Q and the limitation of homogenisation theories to long wavelength 
excitations ||Tj means that most of the interesting effects associated with PCs, which are dynamic in nature and 
therefore correspond to high-frequency bands, are not captured. 


The article pD) described how, by considering perturbations about the standing wave eigenfrequencies of a pe¬ 
riodic structure, a long-scale PDE can be derived that governs the propagation of a scalar modulation function, 
characterising the group behaviour of waves in the structure. Crucially, this approach allows for arbitrarily rapid 
field variation inside the long-scale envelope, which in turn allows us to break free of the quasi-static regime, 
relying only on the conditions of Bloch’s theorem pO) being satisfied. Accordingly, we take advantage of a re¬ 
sult from solid state physics: the dispersive properties of a bulk periodic structure are entirely characterised by 
the frequency dependence of the Bloch wave vector within the irreducible Brillouin zone This result holds 
for general Bravais lattices, and recent work has been done to extend the current theory to such structures | [43) . 
Here, we focus on simple square and cubic lattice arrays for clarity. Since publication of the asymptotic theory in 
2010, further work has been done in which its results have been compared against direct numerical calculations, 
demonstrating that it is capable of capturing dynamic effects unique to high-frequency bands, including all-angle 
negative refraction, dynamic anisotropy and cloaking near Dirac-like cones Q, as well as reproducing the value of 
the reflection coefficient in simple scattering problems |351. This supports the claim that the resulting behaviour 
can be considered that of an effective medium, and justifies our reference to the theory as high-frequency ho¬ 
mogenisation (HEH). Alongside analogous work on discrete lattices pT] , HEH was originally developed for the 
planar Helmholtz equation, which models the frequency-domain response of out-of-plane elastic shear waves, TE 
or TM polarised electromagnetic waves, or pressure waves in a fluid. Eurther work has been done to extend the 
theory to in-plane elasticity, in which the resulting Lame system poses a vector problem in two dimensions 00; 
see also |TT) for further work in elasticity. 

The main result here is to extend the theory to the three-dimensional vector formulation of Maxwell’s equations. 
We shall see that HEH leads to an effective PDE for a scalar modulation function /o, in which the local dispersive 
properties of the structure are captured by a frequency-dependent tensor T. We show that in the low-frequency 
limit, this leads naturally to the expected result of classical homogenisation theory, in which the governing vector 
equation is recovered with effective material parameters. The theory we present is very general in its scope and we 
choose to illustrate its use on two relevant examples, the first being that of guided wave propagation in PCEs, and 


2 










the second being that of dynamic anisotropy upon a structured metafilm. 

The first application we consider is that of dielectric structures that are invariant in one spatial direction. The 
method could also be applied to so-called wire media [38) , but we note that dielectric structures are typical models 
for PCFs and, unlike metallic waveguides, support coupled modes at oblique incidence which are neither polarised 
in the transverse electric (TE) plane nor the transverse magnetic (TM) plane, and therefore demand a full vector 
treatment of Maxwell’s equations. Compared with fully three-dimensional problems, this example has the advan¬ 
tage that the dependence of the field on the third spatial coordinate, say 0 : 3 , can be explicitly factored out as being 
proportional to expli^x^), reducing the problem to a quasi-two-dimensional one (see figure [O) , and therefore is a 
perfect example for testing the theory in a regime where full numerical simulations are tractable. Similar examples 
have been studied by various authors using numerical simulations 1 25][2^ as well as semi-analytic techniques such 
as multipole expansions and our current work complements these, adding further quantitative analysis and 
providing insight via the appealing notion of asymptotic homogenisation. 

The final application we choose is that of guided waves along a structured metafilm of dielectric spheres. Such 
structures |3T), and metasurfaces p 0 |[^ , are of considerable interest as a subclass of electromagnetic metamate¬ 
rials that have advantages in terms of fabrication. We shall consider guided waves within a planar array of spheres 
that decay exponentially normal to the plane and, as such, are described as Rayleigh-Bloch waves 116 51 In 
terms of the three-dimensional Maxwell system, Rayleigh-Bloch waves guided along a linear array of spheres was 
examined in pO) . Due to the computational effort required there have been very limited studies of the waves that 
can exist within planar array of spherical inclusions, with 1 57 1 focussing on scalar acoustics only, and multipole 
methods being used by 


See also 0 ; Rayleigh-Bloch waves are identified but only in limited directions. 
Here we access the full Brillouin zone numerically and identify frequencies for which marked dynamic anisotropy 
occurs, where the effective medium becomes locally hyperbolic, and as a result energy is localised and directed 
along characteristic directions. In simpler two-dimensional situations star-shaped highly directional wave motion 
at specific frequencies for structured media has emerged in experiments and theory in optics 115 and is most 
strikingly seen in discrete mass-spring lattice systems ||^[n 18 371. Here for the full vector Maxwell system the 
three-dimensional simulations are substantial and require many hours of computer time, whereas the asymptotic 
technique identifies these features, provides physical interpretation and finds quantitatively comparable results 
rapidly. 


The article is structured as follows: we begin by setting the scene in terms of the governing equations and necessary 
details regarding periodicity through reciprocal lattices, Brillouin zones and Bloch waves in section]^ Given these 
preliminaries we then move on to developing the asymptotic theory in section]^ and use this as an opportunity 
to clarify some details regarding Dirac-like points (section [TT] l. Furthermore, we confirm that the general theory 
provided does indeed retrieve the quasi-static classical homogenised result (section [T^ . Applications to PCFs 
and then to dynamically anisotropic Rayleigh-Bloch waves on a metafilm are given in sections 00 in which we 
use the asymptotic theory to complement full numerical simulations. Finally, some concluding remarks are drawn 
together in section]^ 


2 Governing equations and reciprocal lattice 

For time harmonic excitations proportional to exp(—iwf), and in the absence of sources. Maxwell’s equations in 
linear media are given by |[^: 


V-D = 0, VxE — jwB = 0, 

V-B = 0, VxH-|- jwD = 0, 


( 2 . 1 ) 
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where E and B are the electric and magnetic fields respectively and are related to the displacement field D and 
magnetizing field H through the constitutive relations 


D = eE, B = nU. (2.2) 

The parameters e = Creo and p = /ir/io are the permittivity and permeability respectively. Our convention is to 
assume that the fields are complex vectors in and the measurable fields are obtained by taking the real part. 
We assume that e and /r are piecewise smooth, spatially varying real parameters that are exactly periodic in each 
direction. Combining the curl equations in with the constitutive relations \2.2\ yields a decoupled equation 
satisfied by the magnetizing field: 


- V X (e-^V X H) + = 0, 


(2.3) 


along with the condition that the component of both E and H parallel to a discontinuity in e or p, must be con¬ 
tinuous. In practice, both e and p are usually piecewise constant, and \23\ has components that are Helmholtz 
equations coupled through boundary conditions between adjacent phases. We note that the magnetic formulation 
( |2.3| l turns out to be more convenient than the symmetric equation for E in the case of dielectrics, because the 
assumption that p is constant ensures that V • H is identically zero, which is useful for the elimination of spuri¬ 
ous modes that may arise in the numerical scheme p5j . Perfect conductors can also be included in the structure, 
imposing a vanishing tangential electric field component at the boundary, corresponding to a vanishing curl of the 
magnetic field parallel to the boundary. This notion is meaningful only in the microwave region of the spectrum, 
which is significantly lower than the plasma frequency of many metals; such cases are of physical interest as this 
is the relevant regime for telecommunications. 


The underlying periodicity of the medium that we seek to homogenise is an important ingredient of the asymptotic 
theory, and has a bearing on the numerical simulations we perform, and some notation is necessary which we 
provide here; this is available in standard texts 112 361 amongst others. Suppose we have a Cartesian co-ordinate 
system with basis vectors Xj for i = 1, 2, 3. We assume that the periodicity of the structure is that of a simple 
cubic array, given explicitly by e(x) = e(x -f 2^mxi + nk -2 -t- pxs]) for m,n,p G Z, and similarly for p(x). 
The primitive lattice vectors 2Z{xi} form an orthonormal set, and the elementary cell C can be chosen as any cube 
of side 21 oriented in accordance with these. Also central to our analysis will be the first Brillouin zone (shown 
in figure [TT| for a simple cubic lattice, along with the analogue for a square lattice in two dimensions), defined in 
terms of the reciprocal lattice vectors ki for i = 1, 2,3, where ki • 2Z{xj} = 27r(5ij JT^ . The relevance of reciprocal 
space here stems from Bloch’s theorem, which states that for an infinite structure with direct lattice vectors {d}, 
the electric field components are quasi-periodic such that H(x-|-d) = H(x) exp(ik d), and similarly for E, where 
in this context the reciprocal vector k is called the Bloch wave vector |Tg. In particular, if the Bloch wave vector 
corresponds to a reciprocal lattice vector, the field has the periodicity of the direct lattice, resulting in a standing 
wave solution. In fact, standing wave solutions are found whenever the Bloch wave vector implies periodicity or 
anti-periodicity in each spatial direction, and it is these points about which the following asymptotic perturbation 
scheme is applied. 


3 The asymptotic procedure 

To adopt a two-scale approach we define the micro-scale position vector ^ = x/l, which is restricted such that 
inside the elementary cell C it takes values G [—1,1] for i = 1, 2, 3. Clearly, variation in e and /i are functions 
of these co-ordinates only. From Bloch theory, we know that the three-dimensional structure admits standing wave 
solutions at discrete eigenfrequencies p6|, corresponding to solutions with Bloch wave vectors at the vertices 
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Figure 2.1: First Brillouin zones for (a) simple square and (b) simple cubic arrays. The triangle in (a), or tetrahedron of (b), whose vertices are labelled are 
the irreducible Brillouin zones, defined as the first Brillouin zones reduced by the point symmetries of the cell, in this case that of the square, cube respectively. 
Standing waves occur at each of the lettered points. 


of the irreducible Brillouin zone. Perturbing the frequency slightly results in long wavelength modulation of the 
Bloch modes, and with this in mind we introduce a macro-scale position vector X = ryx// relative to some fixed 
origin, with 0 < ry <C 1. Provided that the frequency perturbation is not too great, there is a natural separation of 
scales, with X representing the long-scale response of the structure. In the following theory, we are interested in 
the vanishing limit of ry. 



Figure 3.1: Disparate co-ordinate systems in a two-phase triply-periodic structure. The micro-scale co-ordinates (^1,^2,^3) capture the field variation over a 
single elementary cell C — T>i U T>2, whilst the macro-scale co-ordinates {Xi,X2,X^) encode the long-scale behaviour of the medium. 


As is typical in homogenisation theories, the disparity of the length scales associated with ^ and X allows us to 
treat them as independent variables, so the partial derivative operators are expanded using the chain rule as = 
+ v9xi)/l- We label the dimensionless standing wave eigenfrequencies CIq = oJoljc, with c = l/^fioCo the 
speed of light in vacuum, corresponding to solutions satisfying periodic and/or antiperiodic boundary conditions 
on the cell level, explicitly given by 
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(3.1) 


er^V X = ier^V x j = 1,2,3, 

where ± is chosen depending on the vertex of the irreducible Brillouin zone. Following pO) , we seek solutions in 
the form of an asymptotic ansatz using the small parameter rj: 


OO OO 

H = ^77*H„ (3.2) 

z—0 i—0 

which leads to a hierarchy of equations to be solved upwards from the lowest order in rj, with periodic/antiperiodic 
boundary conditions to be applied on the short scale at each level. The leading order system poses an eigenvalue 
problem on the short scale only: 

0(1) : -Vj X er^(|)V 5 x Ho(X, ^) +/r,(|)02Ho(X, = 0, (3.3) 

subject to short-scale periodicity conditions resulting from the two-scale expansion of p.l[ ), and the following 
continuity conditions at any phase interfaces: 


X Hn 


\dT>i 


X n = 0, 


Hn 


dT>i 


X n = 0, 


(3.4) 


where [•] denotes a jump discontinuity and n is normal to the interface dT>i 2 between two adjacent phases. For cer¬ 
tain simple cases, it may be possible to construct analytic or semi-analytic solutions to this problem (for example, 
using multipole expansions for spherical or cylindrical inclusions, analogously to the methods used in |27 40)). 
More often, however, the problem is solved numerically using a finite element package, say pT) , which amounts 
to minimising the residue of the following weak equation: 


'Vj X Ho • Vj X H'* + kVj • MrHo)(V 5 • ^rH'*) - ■ H'*) 


dV = 0. 


(3.5) 


where Hq are the weight functions, <5 is a small positive constant and •* denotes complex conjugation. Although 
for Oq ^ 0 it is usually unnecessary to further impose the divergence-free condition (which follows from taking 
the divergence of p.3| l), we follow the approach of |25| and include the second weak term in (3.5 1 to penalise 
spurious solutions that lie on the lowest branch and do not automatically satisfy this condition. This approach 
alleviates the restriction on geometry or choice of materials. We note that, as we are studying the vanishing limit 
77 , the asymptotic hierarchy is robust even for high-contrast media, a subject that is of great interest in acoustic 
metamaterials, and approached in this context by Zhikov and others (see |56 ^). In practice, a higher contrast 
will limit the frequency range over which we can expect our results to be reliable. Returning to the current problem, 
the general solution to p.3|l is a linear combination of independent vector solutions whose coefficients must be 


allowed to vary on the long scale. Explicitly iFoi(X, = /g’^(X)/ig^''(^, flo) for * = 1, 2j ^nd we sum over 
r = 1, 2, ...,p for an eigenvalue with multiplicity p. Each term is the product of a long-scale scal ar m odulation 
function /g’^^(X) with a vector short-scale that is a known Bloch eigensolution of (3.3i; the aim 




is to find PDEs that the long-scale functions satisfy. Note that the separated-scale component of functions are 
distinguished by lower-case lettering; this convention is followed for the remainder of this paper. We now move to 
the next order in the hierarchy, resulting in an inhomogeneous version of the leading order equation: 


0{r]) : -Vj X x Hi -f /Ti-flgHi = Vj x e,. x Hq + x x Hq - //rfliHo. (3.6) 
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with periodicity conditions as before and continuity conditions given by 


[e/(Vj xHi+VxxHo)]^^^ ^ xn = 0, ^ x n = 0, (3.7) 

at the phase interfaces. Before attempting to solve this, we derive a compatibility condition (see appendix [^, 
whose implication depends on the nature of the eigenvalue in question. The result resembles a p-component Dirac 
equation in long-scale position space: 


(r) 


pnr 

3 


dA_ 

dXi 


+ = 0 , 


(3.8) 


where P"’’ = iUo 


X 


Ar) 


X h, 


(") 


} di/, g”" = [ 

'• 3 Jc 


h (n)* 
0 




Since the operator in (3.31 is Hermitian, we can always choose the eigenvectors hg to be orthonormal such that 


= Snr, and this can be achieved using the Gram-Schmidt procedure. For fixed n, P"" is proportional to the 

j-component of the complex Poynting flux eo X through the cell, which is necessarily zero for standing 
wave solutions satisfying the boundary conditions p.4[ ). This observation is easily proved from the fact that with 
(anti-)periodic boundary conditions it is always possible to normalise the helds such that either or is 
entirely real and the other is imaginary. The value of PJ*’' for n ^ r depends on the nature of the degeneracy un¬ 
der consideration; the hrst and most frequent type encountered is an essential degeneracy, which results from the 
invariance of the system (i.e. the elementary cell and boundary conditions) under a particular symmetry transfor¬ 
mation. For cubic and square arrays, such degeneracies can be shown by symmetry to yield P"’’ = 0. The second 
and more interesting case is that of an accidental degeneracy, which occurs if we continuously tune the parameters 
of the structure until two otherwise distinct eigenvalues occur at the same frequency. Such degeneracies do not 
result from a symmetry group of the cell, and in general P”’’ ^ 0 for n ^ r. Such cases are dealt with in section 
0 


We deduce from (3.8 1 that for an eigenvalue with P”’" = 0 for every n and r then Hi = 0 for non-trivial solutions. 
This results in locally quadratic dispersion bands passing through essentially degenerate eigenvalues (which is not 
necessarily the case for other Bravais lattices, as noted in |fT4)). With this in place, we are in a position to solve 
p.6| l. The solution consists of a complimentary part, which has the same short scale form as the leading order 
solution, and a particular solution modulated by hrst order partial derivatives of /g. The most general form thus 
has components given by -(- /g^^, do)- Exploiting the linearity of 

the system and the independence of the disparate position variables, ( |3.6[ ) then reduces to 3p separate systems of 3 
coupled equations corresponding to each permutation of j, r. These are solved numerically using hnite elements. 
We next move onto the 0{rf‘) system: 


0{rf) : -Vj X Cj. x H 2 -I- = Vj x x x Hi -|- x x Hi 


Cj ^Vx X Vx X Hq — 


(3.9) 


mi 


with periodicity conditions as before and phase boundary conditions: 


[e, xH 2 -fVx xHi)]g^^^ x n = 0, 


= (3.10) 


7 





Here we derive a second compatibility condition (see appendix]^, analogous to that of the previous order. The re¬ 
sult is an effective PDE satisfied by the modulation function /o(X), posed on the long-scale. For each independent 
eigenfunction at Hq, /o is governed by 


T 

1 .' 


dXidX, 


+ njfo = o^T, 


d^fo , - wg) 


dxidx, 


+ 


/o = 0, 


(3.11) 


where the components of the tensor are formed from integrals over the cell involving the leading and first order 
solutions and are given in appendix]^ The equation on the left makes clear the crucial point that the short scale 
is completely absent, having been integrated out and encapsulated by constant tensor whereas the equation 
on the right has been rewritten in terms of dimensional variables, and makes clear that the the long-scale position 
variable X is an artefact of the asymptotic method. In a continuum setting, ( |3.1 l| l permits Bloch wave solutions of 
the form /o = exp(i/« • X/ 77 ) = exp(j/« • x/Z), where the dimensionless quantity k = (k — ko)Z is proportional 
to the difference between the Bloch vector at frequency H and that at Hq, corresponding to the vertex of the 
irreducible Brillouin zone. Substituting this into (|3.1 l|l leads to locally quadratic behaviour of the dispersion bands 


as H 2 = Tij H = Hq -b where for the second relationship we have used the binomial theorem 

to expand the ansatz ( |3.2| [6)) for small k. This serves as a useful verification of the asymptotic method, giving 
expressions for the local behaviour of the dispersion curves, as shown in figure 3.2 for the case of perfectly- 


conducting spherical inclusions. The low-frequency limit, as well as the case of so-called Dirac-like points differ 
from this, both yielding locally linear dispersion and discussed in the following two sections. 



Figure 3.2: (a) Photonic band diagram for a cubic array of perfectly-conducting spheres of radius r — 0.791 in a background of air, as shown in (6). The 
band structure is plotted around the edges of the tetrahedral iri'educible Brillouin zone (c). Solid black curves ai'e from FEM calculation, and red dashed curves are 
asymptotics from equation irm. Note the Dirac-like point at R, which occurs only at this particular value of the radius. 


3.1 Dirac-like points 

The unusual properties assciated with Dirac points in graphene pb) have prompted a surge of interest in the study 
of Dirac-like points in photonics Eig. Such points are characterised by linear crossing of dispersion branches 
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at Brillouin zone vertices, and in terms of the asymptotic theory are associated with non-trivial degeneracies of 
p.3| l. In figurea so-called accidental degeneracy is induced by choosing the radius of the spherical inclusion 
such that two eigenvalues coincide at R point. In contrast, the other degeneracies observed in the band diagram 
are robust under the change of radius or material parameters, and can only be removed by breaking the symmetry 
of the cell; these are called essential degeneracies. 


In the case of an accidental degeneracy (or incidentally for an essential degeneracy in certain non-cubic/square 
lattices), we find that the quantity P"'’ appearing in equation (3.8 1 contains non-zero elements. As a result, the 
dispersion is governed to leading order by (3.8i, which results, unusually, in the appearance of locally linear 
dispersion at non-zero frequencies at a Brilouin zone vertex. It has been noted by various authors gig, that 
unusual wave guiding effects, as well as perfect transmission and cloaking can be observed in some instances 
around Dirac-like points. Such behaviour is considered in detail by Chan et al in p?) , where it is demonstrated 
using an effective medium theory | |59) that under certain conditions the dispersion can be mapped to that of a zero 
refractive index material. This is shown to occur if one can induce an accidental degeneracy between two triply 
degenerate eigenvalues at T point, which is then demonstrated to be possible in a core-shell structure with perfectly 
conducting inclusions. The local dispersion in such cases is characterised by the crossing of two linear bands, with 
relatively flat quadratic bands passing through the middle, as seen in figure [33l 




Figure 3.3: (a) Dirac-like dispersion in a core-shell structure as demonstrated in tlil Here the air-filled elementary cell (fi) contains a spherical shell of 
permittivity 12 surrounding a perfectly conducting core. The effective equation governing the (repeated) linear branches is 0.054/o,Xi + 0.054/o,x2X2 “I" 


0.054/o,x 3X3 + ^i/o = 0, resulting from 


3.12 


In order to extract the relevant asymptotics in the case of a Dirac-like point, we differentiate ( |3.8| l once with respect 
to Xi, and then by self-substitution obtain 


D 


dXidX, 


+ = 0 , 


(3.12) 


where 

J 


and 


= ■ Substituting the Bloch-wave solution (X) = exp(iK • 

X/p) for fixed small k leads to a homogeneous matrix equation that can be solved numerically to give the linear 
asymptotics, as has been done in figure [T2| but in general it is not possible to obtain separate equations governing 
the evolution of the distinct modes. This is consistent with the assertion in that in general no effective medium 
description can be applied satisfactorily at Dirac-like points. Under certain circumstances, however, it does turn 
out to be possible; one such example being that of figure 3.3 with (3.12|i decoupling to four identical isotropic 
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equations; 0.054/o ,XiXi +0.054/o,x2X2 +0-054/o,x3jf3 + ilffo = 0, governing the linear bands, with a final pair 
yielding = 0, and hence with quadratic dispersion governed by (3.11 1 . For the modes with linear dispersion, 
we deduce an effective Helmholtz equation for the Bloch wave, written in dimensional quantities as V^/o + 
n^ffeoHo{uj'^ — Wg)/o = 0, with = eofio{u}'^ — Wg)Z^/0.054 — )• 0 as w — )■ Wq. 


3.2 The low-frequency long-wave limit in dielectric media 


One surprising feature of HFH, when applied to a vector problem, is that the effective PDFs or systems of PDFs 
are scalar equations that do not mirror the form of the governing vector system of ( |2.3[ ). This observation leads 
to an apparent inconsistency with classical homogenisation theory in the low-frequency regime, where the result 
is a vector equation resembling (2.3 i but with effective tensors playing the role of e and fj, |[29). In fact, using the 


current formulation, the effective permittivity and permeability are embedded in the system of equations (B.2i, 
which correctly describes the low-frequency linear asymptotics, but the two tensors cannot in general be extracted 
individually. However, in the commonly-occurring case of purely dielectric media, to which we now specify, the 
effective permeability /r, = 1, and the effective permittivity tensor can then be extracted, as we now proceed to 
show. 

Classical homogenisation theory is quasi-static and limited to dealing with situations where the harmonic frequency 
is sufficiently low that the field on the cell level is characterised by a small perturbation to an otherwise static field. 
This l eads to dispersionless straight lines emanating from the F point at zero frequency such as those seen in 
figure 3.2 Formally, this occurs when = Oijf), corresponding to the substitution Hg = 0 into the leading 
order system, along with the divergence-free condition Vj • Hg = 0, which is no longer automatically satisfied 
by solutions of p.3| l. Following the HFH procedure of section]^ we deduce that the leading order eigensolution 
is a constant vector field, and Hg = 0 is thus to be treated as a repeated eigenvalue with multiplicity p = 3, 
corresponding to the three spatial directions. After some algebra, the resulting system of coupled effective PDFs 
can indeed be written as a single vector equation, which is the familiar result from classical homogenisation theory 
given in the context of dielectric media in | 


- V X X Ho) + w^eoPoHo = 0, (3.13) 

where and Hg is the leading order magnetic field. Here (•) denotes a cell-averaged 

quantity, and the correction tensor ^ has components depending on the geometry of the cell. To prove this, we 
note that in the quasi-static limit Hg = 0, the leading order solutions can be chosen without loss of generality 



rpnr 



SjnSir + d^h^inj “ ^nh^iij) 


(3.14) 


It is straightforward to show that the first two terms in the integrand satisfy the necessary symmetries, which are 
inherited directly from the —Vx x Vx x (•) operator in (B.l i. It is also trivial to note that the remaining two 
terms are antisymmetric in n and i. The only non-trivial symmetry that remains is antisymmetry of the second two 
terms with respect to j and r. In order to prove this, we consider the inhomogeneous problem (|3.6|l: substituting 
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(r) (r) 

in the constant leading order solutions, along with the candidate particular solution Hu = /q x ^uj to the 
following coupled problems for i,j,r G {1,2, 3}: 


dk{e, - Okie, ^ - 5irdje^ \ 


lij' ‘-'KV^r Ikj 

subject to boundary conditions at any phase interface given by 


(3.15) 


.- 1 ; 


dkh^ujUk - t, '^dih'[Xnk 




ikj' 



-1 

— 


' 2 ^ 1,2 



•Di 


{SijTlj- SirTlj')^ 


(3.16) 


where is the component of the unit normal n in the Xi-direction and [•] denotes a jump discontinuity. The 
two-scale expansion of the divergence-free condition Vj • Hi = —Vx • Hq must also be imposed, which gives 


f{r) Ar) 


= -/, 


(r) 


Q X ^jr- In the case that j ^ r, this condition is homogeneous, and because the right hand sides 


. . J I (r) (7") 

of both (|3.15|l and (|3.16| are antisymmetric in j, r, we deduce that On the other hand, for j = r 


the divergence-free condition yields h\^- ^ = —1. Integrating this over the cell C and applying the divergence 
theorem and periodicity leads to a contradiction “0 — 1” that can only be avoided if /q x ~ hence 


derived the long-scale divergence-free condition Vx • Hq = 0 , and with this in place we can show for j = r 
that h^^lj = 0. In doing so, T is 
effective permittivity as follows: 


(r) 

that h\ -J = 0. In doing so, T is deduced to contain the correct symmetries, and we derive an expression for the 


rpnr _ -l,hom 

-Lij — ^mk^jrltf^i 


— l,hom 
^kl 


{2\y 


I ^ni kSjrlTij 


inr 

ij ■ 


(3.17) 


For each permutation of k, I, we only need evaluate one component of T. It is natural to split the resulting inverse 
permittivity tensor into a sum of two parts, the hrst resulting from the Kronecker delta terms in the integrand of 
)5ij, and the second from the two remaining terms, which we refer to as the correction tensor 


Tyj', given by ^ 

e“7. Putting all this together we recover (|3.13| and hence the quasi-static homogenisation of Maxwell’s equations. 


4 Application to PCFs 


We now turn our attention to the case of dielectric cylinders of inhnite length and constant permittivity Cc aligned 
in the xa-direction and embedded in a matrix phase with permittivity e^n- This model is appropriate for both holey- 
type PCFs, in which the cylinders are air holes in a dielectric background, as well as ARROW-type PCFs, in which 
the cylinders have a higher refractive index than the background. As discussed in section [T] the a; 3 -dependence of 
the solution can be factored out as oc exp(i/3a;3) where /? is the propagation constant of the radiation. 


For a fixed value of (3, the resulting Bloch wave structure is projected onto the (xi , X 2 ) plane and the effective PDFs 
will govern propagation in these directions. Tuning /?, along with the geometric and material properties of the cell, 
can then be used to induce particular features, such as partial stop-bands or Dirac-like points in the transverse 
plane. In hgurejA^ we show the band diagrams for one such conhguration of air holes in a dielectric background 
at two different values of /3 chosen such that a transverse stop-band and a Dirac-like point are exhibited. 


In order to obtain the quasi-planar analogues of (3.8 1 and ( |3.1 l| l, we simply make the substitutions d/dX^ -G 0 
and d/d ^3 —?> if3, and the corresponding integrals become surface integrals over the two-dimensional cell. It is 
also straightforward to check for /3 = 0 that our theory reduces to that of planar HFH, as published in j 
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Figure 4.1: Reduction of a typical fibre-like structure to two dimensions. As shown in (a), dielectric cylinders of permittivity tc are embedded in a homo¬ 
geneous matrix phase with permittivity The resulting problem is projected onto the (xi, X 2 ) plane, and the angle of propagation depends on the ratio of the 
parameter 0 and the planar Bloch wave vector, which characterises the phase shift across an irreducible cell (the dashed square in (b)) 



Figure 4.2: Transverse photonic band diagrams for the lowest 5 bands for a square array of cylindrical air holes of radius r = 0.75/, where 21 is the pitch 
of the an'ay, in a dielectric background with relative permittivity Cr = 6. In (a), the parameter/3 = 0.89/c/has been tuned to induce an accidental degeneracy, 
resulting in a Dirac-like point at M. In (b), 0 has been increased to 3/1, causing the higher-frequency bands to separate from the lowest two, so the PCF exhibits 
a partial band gap in the transverse plane. Note that different scales have been used on the vertical axes. 
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4.1 Transverse stop bands and localised defect modes 


So far we have focused on the application of HFH to propagating ‘slow-modes’ inside a bulk periodic structure, 
and the corresponding modes in the quasi-planar case propagate in the transverse plane as well as in the axial 
direction. Alternatively, we may consider a structure exhibiting a transverse stop band at a particular value of 
p. If we introduce a finite defect of some kind into the structure, it may be possible to excite a mode whose 
frequency lies within the stop band, and hence decays evanescently in the surrounding medium. If the frequency 
lies in the vicinity of one of the eigenfrequencies flo, the exponentially decay is then governed by ( 3.11| l. In figure 
|4.3| we show one such mode, where a defect is induced by removing one air cylinder from a large array of cells 
corresponding to those of fig |4.2{ 6). Such a model is highly relevant for real-life applications of PCFs, and one can 
envisage the extension of such work giving rise to the notion of ‘effective fibres’. 



Figure 4.3: Magnetic field norm for a localised defect mode at = 1.690, induced by removing one air cylinder from a large array of cells (of which a portion 
are shown here), corresponding to those of fig |4.2| h). The dominant effective equation is 1.01 + 0.098/o,X2X2 +^ 2/0 — 0, where Oq — 1.707, 

corresponding to the third eigenvalue at X, which leads to directed leakage in the xi-direction, and hence also in the ai 2 -direction due to the symmetry point at the 
opposite vertex of the first Brillouin zone. The decay of the mode in these directions is then very well approximated by an exponential form exp( —aai) (the red 
dashed line in (b)), where the decay rate a = 0.24/1 follows from the PDE. 


5 Application to Rayleigh-Bloch modes 

For various wave systems, it is well known that periodic arrays of defects in otherwise homogeneous media can 
support localised modes. Such modes exhibit Bloch wave quasi-periodicity in the array, and decay evanescently in 
the surrounding medium, reminiscent of Rayleigh waves at the surface of solids. It is therefore natural to describe 
them as Rayleigh-Bloch modes, and their existence is ubiquitous in systems where the dimensionality of the array 
is at least one less than that of the wave system in question. Rayleigh-Bloch waves are distinguished from “pure” 
surface waves (e.g. Rayleigh waves. Lamb waves, etc.) in that they can propagate along surfaces that exhibit some 
material or geometric periodicity and which, in the absence of this periodicity, do not support surface waves. In 
this sense they are similar to spoof surface plasmons ^9\ . In two- and three-dimensional systems, linear arrays 
may be employed as diffraction gratings, and have been studied in the context of water waves and acoustics as 
solutions to the Helmholtz equation | |^[5^|57) , surface plasmons | |42||45) , as well as coupled elastic waves HD- 
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In three-dimensional systems, such as the vector Maxwell system, Rayleigh-Bloch modes can exist for linear or 
planar arrays 1 40 1 , the latter allowing for the possibility of a range of in-plane effects to be observed. 



Figure 5.1: A Doubly-periodic planar square array of dielectric spheres. For the purpose of computation, the spheres have radii of 0.8Z, where the pitch of the 
array is 21. The spheres are embedded in free space and have a relativity permittivity Cr = 20. 


The asymptotic technique we have developed is particularly well-suited to studying electromagnetic Rayleigh- 
Bloch systems, as the corresponding full numerical computations can be prohibitively demanding, particularly 
for fully coupled three-dimensional problems of the type examined here. In order to assess and demonstrate the 
efficacy of the asymptotic technique developed in the preceding sections, we consider an infinite doubly-periodic 
planar array of dielectric spheres embedded in free space, as depicted in figure [5T| Once again, the theory is largely 
unaltered from the three-dimensional case; we simply make the substitution d/dX^ —0 in (3.8 1 and ( |3.11[ ), and 
extend the cell in the 0:3 direction as shown in figure |5^6). The corresponding photonic band diagram is shown 
inl^a). 


5.1 Planar dynamic anisotropy & hyperbolic metafilms 

The band diagram of figure [5^ and has several interesting features including a planar band-gap as well as flat bands 
associated with so-called slow waves. A particularly fascinating feature is the change in curvature of the lowest 
branch at lattice point X in the reciprocal space, shown in the inset of figure 5.2 a) along with the corresponding 
local asymptotics. Such saddle points are associated with dynamic anisotropy l^[Tg| 8 ) and hyperbolic materi¬ 


als 15011541 ; this type of media, when excited by an appropriate multipole source, direct the electromagnetic energy 
along clearly defined characteristic lines. The existence of these saddle points for the system under consideration 
suggests the interesting possibility of dynamic anisotropy in a Rayleigh-Bloch wave setting, creating an effective 
hyperbolic metasurface. 

Using the asymptotic theory developed earlier, we first obtain the leading order solution for the field when the array 
is excited by a dipole source at a frequency close to that of the saddle point on the dispersion surface (Oq = 0.8071). 
This involves solving the leading and first order cell problems as outlined in section]^ and the solvability condition 
for the second-order problem then generates a PDE, of the same form as (|3.11|l, for the long-scale field; 


rj. d^fo 
dxidx^ 




/o = 0, 


(5.1) 


where T is now a rank-2 tensor. In this case, we find that T = diag(—0.19,0.11) is diagonal and T 11 T 22 < 0, 
yielding a hyperbolic PDE on the long-scale. The leading order asymptotic field for the forced problem is then a 
superposition of the solution to this PDE, whose forcing must take into account the symmetry of the cell problem, 
along with its symmetrical counterpart excited at point Y in the Brillouin zone. The result is most striking when 
the electric field is plotted, as in figure [53a where since the hyperbolic mode is primarily out-of-plane electric. 
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Figure 5.2: (a) In-plane photonic band diagram for a square array of high-permittivity (cr — 20) dielectric spheres in vacuum, as shown in figurejsjJ The 
radius of the spheres is 0.8/, where 21 is the pitch of the array, corresponding to the elementary cell (b). Here the cell is extended in the ^ 3 -direction, and perfectly 
conducting boundary conditions are applied at its top and bottom of the cell. The inset of (a) shows an enlarged portion of the diagram in which the local behaviour 
is governed by the hyperbolic effective equation —0.19/o,XiXi + 0 . 11 /o,X 2 X 2 d" ^ 2/0 ~ leading to the asymptotics shown by red dots. 


hence we only include the i? 3 -component of the field. Here, the non-dimensional angular frequency of excitation 
is O = 0.80705, which lies close to the resonant frequency of Oq = 0.8071. Both the cell and long-scale problems 
were solved using the commercial finite element package Comsol Multiphysics®. 


For comparison the above problem is also treated, purely numerically, using finite elements. In particular, a 37 x 37 
planar array of spheres is considered; the finite cluster is surrounded by regions of perfectly matched layers in order 
to simulate an infinite domain. Moreover, the computational cost of the full finite element simulation is lessened 
by making full use of the available symmetries, thereby reducing the computational domain by a factor of 8. The 
array is then excited, at the centre of the array, by a dipole source of unitary magnitude and oriented along the 
a; 3 -axis. The component of the electric field from the full numerical simulation is shown in figure 5.3b and 
should be compared with the asymptotic field shown in figure 5.3a For comparison, we also plot the field along 
the line (x, y, z) = {x, x, 0) and passing through the dipole source (see figure [T^ . 

Finally, figure [S^ shows the component of the electric field both in the plane of the array and, over two slices, 
perpendicular to the array. The figure illustrates the decay of the field in the direction perpendicular to the array in 
addition to the dynamic anisotropy exhibited in the plane of the array. The novel exhibition of dynamic anisotropy 
on a structured metafilm is an exciting effect with several potential applications in the guiding of surface waves; this 
effect allows not only the confinement of waves to a surface or interface, but also the control of their propagation 
within the metafilm itself. These effects occur, by necessity, at frequencies where the wavelength is comparable to 
the size of the microstructure where traditional (long-wave) homogenisation theories are no longer valid. However, 
as demonstrated by figure [53j the asymptotic homogenisation scheme developed here is capable of accurately and 
conveniently describing such effects. 


Although we have not carried out a comprehensive analysis of the computational cost of the asymptotic method, 
it is illuminating to consider the following: The full finite element simulation, using all the available symmetries. 
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(a) Asymptotic analysis (b) Finite element simulation 



(c) 


Figure 5.3: The z-component of the electric field for the planar array obtained from (a) the asymptotic analysis and (b) the full finite element simulations. In 
both cases the colour scale is linear running for minimum (blue), through zero (green) to maximal (red). Pail (c) shows the 2 -component of the electric field plotted 
along a line passing through the points x — (0,0,0) and (1, 1,0). The solid curve represents the field from the full finite element simulation, whereas the dashed 
curve coiTesponds to the leading order solution from the asymptotic analysis. The non-dimensional angular frequency of excitation is — 0.80705, which is 
close to the resonant frequency of f2o = 0.8071. 
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Figure 5.4: Slices through the array showing the S 3 component of the electric field and illustrating the decay of the surface wave in the direction perpendicular 
to the array together with the in-plane dynamic anisotropy. Once again, the colour scale is linear running for minimum (blue), through zero (green) to maximal 
(red). 


required approximately 16.8 million degrees of freedom and took 11.5 hours to solve on a dedicated machine using 
77GB of RAM and 16 processor cores running at 2.5GHz. In comparison, the numerical element of the asymptotic 
scheme used around 830 thousand degrees of freedom, 3.7GB of RAM and required 2.7 minutes to solve on a 
laptop with 2 cores running at 2.5GHz. 


6 Concluding remarks 

We have demonstrated that in the vicinity of Brillouin zone vertices for a periodic structure, photonic modes gov¬ 
erned by the vector Maxwell system give rise to second-order scalar PDEs, whose solutions govern the propagation, 
or decay, of waves inside the structure. In some cases, these PDEs can be used straightforwardly to predict quali¬ 
tatively and quantitatively the long-scale behaviour of a periodic photonic structure, which has been demonstrated 
in the context of a forced problem in a planar PC, as well as an eigenvalue problem for a localised mode in a PCE. 

A particularly interesting physical example was considered in section wherein the asymptotic homogenisation 
theory was applied to a fully coupled three-dimensional problem for the Maxwell system. Using a planar array 
of dielectric spheres we were able to create an effective hyperbolic metafilm which supports surface waves that 
also exhibit dynamic anisotropy over the surface. One can envisage many applications where the ability to guide 
electromagnetic waves along a surface whilst also controlling the propagation over the surface itself would be 
useful. 

A natural extension of this work is to consider more general Bravais geometries that are not necessarily square/cubic. 
In fact, such work has recently been completed for the simpler case of TE-polarised radiation | |43) , and it is clear 
that the majority of the theory presented here carries through mostly unaltered. A second extension would be to 
consider structures with a wider range of material properties, for example those experiencing dielectric loss, or 
perhaps more interestingly real metals, whose frequency-dependent permittivity can be described by the Drude 
model. 
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A The 0 { t ]) compatibility condition 


For each n S {1, 2, ...p}, we first take the scalar product of (3.6 1 with and subtract from the scalar product 


of (3.3 I* with Hi//q”\ The resulting equation reads 


• (Vj X X Hi) - Hi • (Vj x x h),”)*) 


= h 


(n)* 


( - Vf X e ^ Vx X Ho - e, ^ Vx X X Ho + Pr^lUo ). 


(A.l) 


Using the vector identity A • (V x V x B) — B • (V x V x A) = V • (B x V x A — A x V x B) |24 , the left 
hand side of this is 


• (Hi X Vj X X Vf X Hi) + (V^e-^) • (Hi x x x x Hi), (A.2) 

which is simply the divergence of the quantity (Hi X Vj X X X Hi). Integrating this over 

the elementary cell C and applying the divergence theorem leaves us with an integral over the surface of the cell 

fnl 

which vanishes by periodicity/antiperiodicity of Hi and , along with a non-vanishing contribution at the phase 

boundaries given by 


f [e.-i{h(")*-(nx Vj xHi)-Hi.(nxV5 xh),")*)}]d5, (A.3) 

J UOX^l ,2 


where again square brackets denote a jump discontinuity. The boundary conditions in (3.41 and ( |3.7| l ensure that 
the only non-zero contribution comes from the first term above, and we are left with the following equation; 


lud-Di 




(nx Vj X Hi)]dS' 


= / • { - Vj X X Ho - er'Vx x x Hq + p,fI?Ho}d'P, 


(A.4) 


for n = 1, 2, ...,p. Using the cyclic property of the scalar triple product, the first term of the integrand on the right 
is Vj • (h),")* X X Ho) - x h[,">*) • (Vx x Ho). Applying the divergence theorem to the first 

resulting term yields a surface integral which, with the help of the boundary conditions in p.4| i and p.7[ ), exactly 
cancels the surface integral on the left hand side of (|A.4|i, leaving the equation 


0 = ^ - er' (V^ X h^”^*) • (Vx X Ho) - • Vx x x Ho + • Ho}dU. (A.5) 

Now, substituting the solution Ho = and noting that x'Ve X = —iUoeo*'"^ leads to the system of 

equations (|3.8|). 


18 






B The 0{rf‘) compatibility condition: effective PDEs 


We now derive a compatibility condit ion f rom p.9[ ) following a similar methodology to th at of appendixwe 
begin by taking the scalar product of (3.9 1 with subtract from the scalar product of (3.3 1 * with H 2 //q"^ 

and integrate the resulting equation over the elementary cell C. The left hand side vanishes like the first order 
equivalent, leaving 


0 = ^ - eri(V 5 X • (Vx X Hi) - • Vx x Vj x Hi 




- e. h, 


0 


Vx X Vx X Ho + • HojdV. 


(B.l) 


Using the cyclic property of the scalar triple product, it is straightforward to show that the contribution of the 
complementary part of Hi in the above equation is zero. Changing to tensor notation, and making use of the 
contracted epsilon identity eijkSkim = SuSjm — SimSji, we are led to the following system of effective PDEs; 


rj^l 




-■) 


dX,dX, 


+ = 0 


rpnr 


dXSX, 


+ = 0 , 


where , with Q'*’' = 


(n)H 

Ofc 


h^^k^V and 


(B.2) 


rpnr 


-1 J r 

Ef \ ^ij\k ^Ok ^Oj ^Oi 


+h^.dkh^r-h^^Ah^r+h, 






,(")♦ , X'^)* 


^Ok 


''Ofc 




■dU. 


(B.3) 


In the case of a distinct eigenvalue, the above system reduces to a single equation with n = r = 1. For a repeated 
eigenvalue, the substitution /g’^^(X) = exp(j«; • X/yy) yields a homogeneous matrix equation whose solution 
gives the quadratic asymptotics, but in fact the independent Bloch modes are non-interacting in the sense that they 
can be decoupled as: 


= 0, (B.4) 

where M is matrix of eigenvectors shared by the the matrices r^, and hence are all diagonal matrices. Pre¬ 
multiplying by M~^, we are then left with the decoupled system 






(r) 


dx,dx, 


+ = 0 , 


(B.5) 


where = (M We can write the general leading order solution as a linear combination of these 


10 

decoupled modes Hq = with hg"'' = M’'”hg ^ Note that the tildes have been omitted in equation (3.111 

as it is assumed that the decoupling is done automatically. 






19 









References 


[1] A. Alu. First-principles homogenization theory for periodic metamaterials. Phys. Rev. B, 84:075153, 2011. 

[2] T. Antonakakis, R. V. Craster, and S. Guenneau. Asymptotics for metamaterials and photonic crystals. Proc. 
R. Soc. A, 496, 2013. 

[3] T. Antonakakis, R. V. Craster, and S. Guenneau. High-frequency homogenization of zero frequency stop 
band photonic and phononic crystals. New J. Phys., 15, 2013. 

[4] T. Antonakakis, R. V. Craster, and S. Guenneau. Homogenization for elastic photonic crystals and dynamic 
anisotropy. J. Mech. Phys. Solids, 71:84-96, 2014. 

[5] M. V. Ayzenberg-Stepanenko and L. I. Slepyan. Resonant-frequency primitive waveforms and star waves in 
lattices. J. Sound Vib., 313:812-821, 2008. 

[6] N.S. Bakhvalov and G. Panasenko. Homogenisation: averaging processes in periodic media. Springer- 
Verlag, Diisseldorf, 1989. 

[7] P. A. Belov and C. R. Simovski. Homogenization of electromagnetic crystals formed by uniaxial resonant 
scatterers. Phys. Rev. E, 72:026615, 2005. 

[8] H. Benisty. Photonic crystals: New designs to confine light. Nat Phys, 1(1):9-10, 10 2005. 

[9] A. Bensoussan, J. Lions, and G. Papanicolaou. Asymptotic analysis for periodic structures. North-Holland, 
Amsterdam, 1978. 

[10] F. Bloch. Uber die Quantenmechanik der Elektronen in Kristallgittern. Z. Physik, 52:555-600, 1928. 

[11] C Boutin, A Rallu, and S Hans. Large scale modulation of high frequency waves in periodic elastic compos¬ 
ites. J. Mech. Phys. Solids, 70:362-381, 2014. 

[12] L. Brillouin. Les electrons dans les metaux et le classement des ondes de de Broglie correspondantes. 
Comptes Rendus Hebdomadaires des Seances de lAcademie des Sciences, 191:292, 1930. 

[13] L. Brillouin. Wave Propagation in Periodic Structures. Dover, 1946. 

[14] C. T. Chan, F. Huang, F. Liu, and Z. H. Hang. Dirac dispersion and zero-index in two dimensional and three 
dimensional photon and phononic systems. Prog. Electromagn. Res. B., 44:163-190, 2012. 

[15] D. N. Chigrin, S. Enoch, C. M. Sotomayer Torres, and G. Tayeb. Self-guiding in two-dimensional photonic 
crystals. Optics Express, 11:1203-1211,2003. 

[16] D. J. Colquitt, R. V. Craster, T. Antonakakis, and S. Guenneau. Rayleigh-Bloch waves along elastic diffrac¬ 
tion gratings. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sci¬ 
ences, 471(2173), 2014. 

[17] D. J. Colquitt, 1. S. Jones, N. V. Movchan, and A. B. Movchan. Dynamic anisotropy and localization in elastic 
lattice systems. Waves in Random and Complex Media, 22(2): 1-17, December 2011. 

[18] R. V. Craster, T. Antonakakis, M. Makwana, and S. Guenneau. Dangers of using the edges of the Brillouin 
zone. Physical Review B, 86:115-130, 2012. 

[19] R. V. Craster, T. Antonakakis, M. Makwana, and S. Guenneau. Dangers of using the edges of the Brillouin 
zone. Physical Review B, 86(11): 1-6, September 2012. 


20 



[20] R. V. Craster, J. Kaplunov, and A. V. Pichugin. High-frequency homogenization for periodic media. Proc. R. 
Soc. Lond, 466:2341-2362, 2010. 

[21] R. V. Craster, J. Kaplunov, and J. Postnova. High-frequency asymptotics, homogenisation and localisation 
for lattices. Q J Mechanics Appl Math, 63:497-519, 2010. 

[22] S.A. Cummer and D. Shurig. One path to acoustic cloaking. New J. Phys., 9:45, 2007. 

[23] G. Dupont, M. Farhat, A. Diatta, S. Guenneau, and S. Enoch. Numerical analysis of three-dimensional 
acoustic cloaks and carpets. Wave Motion, 48:483^96, 2011. 

[24] M. Femandez-Guasti. Green’s second identity for vector fields. ISRNMath. Phys., 2012:973968, 2012. 

[25] S. Guenneau, A. Nicolet, and F. Zolla. A finite element formulation for spectral problems in optical fibers. 
COMPEL, 45(1):120-131, 2001. 

[26] S. Guenneau, C.G. Poulton, and A.B. Movchan. Oblique propagation of electromagnetic and elastic waves 
for an array of cylindrical fibres. IEEE Trans. Magn., 38(2):1261-1264, 2002. 

[27] S. Guenneau, C.G. Poulton, and A.B. Movchan. Oblique propagation of electromagnetic and elastic waves 
for an array of cylindrical fibres. Proc. R. Soc. Lond. A, 459:2215-2263, 2003. 

[28] S. Guenneau and F. Zolla. Homogenization of three-dimensional finite photonic crystals. Prog. In Elect. Res., 
27:91-127, 2000. 

[29] S. Guenneau, F. Zolla, and A. Nicolet. Homogenization of three-dimensional photonic crystals with hetero¬ 
geneous permittivity and permeability’. Waves in Random and Complex Media, 17(4):653-697, 2007. 

[30] C. F. Holloway, E. F. Kuester, J. A. Gordon, J.O’Hara, J. Booth, and D. R. Smith. An overview of the 
theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials. IEEE Antennas 
Propagat. Mag., 54:10-35, 2012. 

[31] C. F. Holloway, D. C. Fove, E. F. Kuester, J. A. Gordon, and D. A. Hill. Use of generalized sheet transition 
conditions to model guided waves on metasurfaces/metafilms. IEEE Trans. Antennas Propagat, 60:5173- 
5186, 2012. 

[32] X. Huang, Y. Fai, Z. Hong Hang, H. Zheng, and C. T. Chan. Dirac cones induced by accidental degeneracy 
in photonic crystals and zero-refractive-index materials. Nat Mater., 2011. 

[33] J.D. Jackson. Classical electrodynamics. Wiley, 1962. 

[34] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade. Photonic Crystals, Molding the Flow of 
Light. Princeton University Press, Princeton, second edition, 2008. 

[35] F. M. Joseph and R. V. Craster. Reflection from a semi-infinite stack of layers using homogenization. Wave 
Motion, 45:145-156, 2015. 

[36] C. Kittel. Introduction to solid state physics. Wiley, New York, 1996. 

[37] R. S. Fangley. The response of two-dimensional periodic structures to point harmonic forcing. J. Sound Vib., 
197:447^69, 1997. 

[38] F. Femoult, M. Fink, and G. Ferosey. Revisiting the wire medium: an ideal resonant metalens. Waves in 
Random and Complex Media, 21:591-613, 2011. 


21 



[39] C M Linton, R Porter, and I Thompson. Scattering by a semi-infinite periodic array and the excitation of 
surface waves. SIAM Journal on Applied Mathematics, 67(5): 1233-1258, 2007. 

[40] C. M. Linton, V. Zalipaev, and 1. Thompson. Electromagnetic guided waves on linear arrays of spheres. Wave 
Motion, 50:29^0, 2013. 

[41] COMSOLltd. Comsol multiphysics 5.0, 2014. 

[42] Stefan a. Maier, Steve R. Andrews, L. Martm-Moreno, and F. J. Garcia-Vidal. Terahertz surface plasmon- 
polariton propagation and focusing on periodically corrugated metal wires. Physical Review Letters, 
97(17):!^, 2006. 

[43] M. Makwana, A. Antonakakis, B. Maling, S. Guenneau, and R.V. Craster. Wave mechanics in media pinned 
at Bravais lattice points. arXiv: 1505.02889, 2014. 

[44] A. A. Maradudin, editor. Structured Surfaces as Optical Metamaterials. Cambridge University Press, 2011. 

[45] M. Navarro-Cia, M. Beruete, S. Agrafiotis, F. Falcone, M. Sorolla, and S. Maier. Broadband spoof plas- 
mons and sub wavelength electromagnetic energy confinement on ultrathin metafilms. Optics express, 
17(20):18184-18195, 2009. 

[46] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim. The electronic properties of 
graphene. Rev. Modern Phys., pages 109-162, 2009. 

[47] O.A. Oleinik, A.S. Shamaev, and G.A. Yosifian. Mathematical problems in elasticity and homogenization. 
North-Holland, Amsterdam, 1992. 

[48] G. G. Osharovich, M. V. Ayzenberg-Stepanenko, and O. Tsareva. Wave propagation in elastic lattices sub¬ 
jected to a local harmonic loading. 11. Two-dimensional problems. Continuum Mechanics and Thermody¬ 
namics, 22(6-8):599-616, August 2010. 

[49] J. B. Pendry, L. Martm-Moreno, and F. J. Garcia-Vidal. Mimicking surface plasmons with structured surfaces. 
Science (New York, N.Y.), 305(5685):847-848, 2004. 

[50] A. Poddubny, 1. lorsh, P. Belov, and Y. Kivshar. Hyperbolic metamaterials. Nature Photonics, 7(12):948-957, 
2013. 

[51] R. Porter and D. V. Evans. Rayleigh-Bloch surface waves along periodic gratings and their connection with 
trapped modes in waveguides. J. FluidMech., 386:233-258, 1999. 

[52] R Porter and D. V. Evans. Embedded RayleighBloch surface waves along periodic rectangular arrays. Wave 
Motion, 43(l):29-50, November 2005. 

[53] P. St. J. Russell. Photonic crystal fibers. Science, 299:358362, 2003. 

[54] P. Shekhar, J. Atkinson, and Z. Jacob. Hyperbolic metamaterials: fundamentals and applications. Nano 
Convergence, 1(1):14, 2014. 

[55] R. A. Shore and A. D. Yaghjian. Traveling waves on two- and three-dimensional periodic arrays of lossless 
scatterers. Radio Sci., 42:RS6S21, 2007. 

[56] J.D. Smith. Application of the method of asymptotic homogenization to an acoustic metamaterial. Proc. R. 
Soc. Lond, 467:3318-3331, 2011. 


22 



[57] I. Thompson and C. M. Linton. Guided surface waves on one- and two- dimensional arrays of spheres. SIAM 
Journal on Applied Mathematics, 70(8):2975-2995, 2010. 

[58] C. H. Wilcox. Scattering theory for diffraction gratings. Springer-Verlag, 1984. 

[59] Y. Wu, J. Li, Q. Zhang, and C. T. Chan. Effective medium theory for magnetodielectric composites: Beyond 
the long-wavelength limit. Phys. Rev. B., 74(085111), 2006. 

[60] V.V. Zhikov. On an extension of the method of two-scale convergence and its applications. Sb. Math., 
191:973-1014, 2000. 

[61] V.V. Zhikov, S.M. Koslov, and O.A. Oleinik. Homogenization of differential operators and integral function¬ 
als. Springer-Verlag, Diisseldorf, 1994. 

[62] F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau, and D. Felbacq. Foundations of photonic 
crystal fibres. Imperial College Press, London, 2005. 


23 



